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ABSTRACT 


The NPSATI satellite uses an active magnetic torque rod system, with a 
magnetometer for attitude determination, to maintain 3-axis stabilization, with a slightly 
gravity gradient friendly structure. 

This thesis will examine the performance of three combinations of programs and 
simulation models for the NPSAT1 satellite attitude control system. The models include 
a magnetic control law with a reduced order estimator to generate torque commands to 
achieve spacecraft nadir pointing and a magnetic rate (Bdot) control law to reduce 
spacecraft angular rates. The performances of two Bdot mode switching designs are 
compared. Also, a case is made for the benefits of priming the system’s reduced 


estimator prior to mode switching. 


All of the control methods analyzed appear to be valid control methods to achieve 
three-axis attitude stabilization using only magnetic torquers for active control. The most 
efficient control method analyzed incorporates a hand-off method from a magnetic rate 
(Bdot) control loop to a magnetic control loop. The results of this analysis indicates that 
the best use of this method is to perform the Bdot hand-off following the achievement of 


a predetermined combined angular rate. 
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EXECUTIVE SUMMARY 


The development of the attitude control system (ACS) for the NPSAT1 satellite 
was driven primarily by one object, namely, student learning through analysis and design 
of an actual spacecraft. Another concern driving the NPSAT1 satellite’s development 
effort was the motivation to produce a small inexpensive research satellite composed of 
COTS (Commercial Off The Shelf) components and space flight experiments. This 


naturally led toward a system that was sensitive to the anticipated low power budget. 


During the effort to keep the above-mentioned concerns in the forefront of the 
attitude control system’s design process, the design team elected a control system that 
was to maintain three-axis stabilization using only magnetic torquers and magnetometers. 
Barry Leonard (Naval Postgraduate School Adjunct Professor) produced two potential 
system models. Both models use a magnetic control law and a reduced estimator to 
determine the spacecraft states and achieve three-axis stabilized nadir pointing using 
magnetic torque rods and a magnetometer. The second model also uses an added 
magnetic rate (Bdot) control law as suggested by the Naval Research Laboratory for 
initial de-tumbling of the spacecraft after deployment. [Ref. 1] 


This thesis examines the performance of the two systems as stand-alone models 
and when used in conjunction with the Bdot control program. One application of the 
Bdot control loop hands off attitude control after a predetermined time interval and a 
second application after a predetermined combined angular rate has been achieved. After 
assigning initial conditions and establishing a test strategy, the models and programs 
were modified to include data capture, storage and plotting algorithms. Due to the need 
to find the lowest ACS power budget possible, in order to allow greater power 
availability for research payloads, this analysis is sensitive to the trade-offs between 
system power consumption and time to achieve three-axis attitude control after launch 
and deployment or after a loss of attitude control due to disturbance forces. The ACS 
needs to achieve and maintain stability within operational requirements (currently nadir 
pointing +/- 10 degrees), yet still maintain the robustness and capability to recover from a 


loss of attitude. The results from each of three combinations of control techniques is 


XV 


summarized and plotted. Nine sets of initial conditions were tested with all three model 


and program combinations. 


All of the control methods analyzed appear to be valid control methods to achieve 
three-axis attitude stabilization using only magnetic torquers and magnetometers. The 
most efficient control method analyzed appeared to be the model that incorporates a 
hand-off method from a Bdot control loop to an estimator control loop. The analysis 
results indicate that the best use of this method is to perform the Bdot hand-off following 


the achievement of a predetermined combined angular rate. 


XV1 


I. INTRODUCTION 


The attitude of a satellite can be controlled by the interaction between the earth’s 
magnetic field and a magnetic moment produced within the spacecraft. Many attitude 
control systems have been designed using magnetic torquers to supplement other active 


or passive attitude control methods. 


The purpose of the NPSAT1 ACS is to prove a satellite’s attitude can be 
maintained solely with magnetic torquers and a 3-axis magnetometer. Magnetic torquers 
rarely fail on orbit and multitudes of satellites have magnetic torque rods supplementing 
their control systems. The attitude control capability of the magnetic torquers, however, 
has limited effectiveness (+/- several degrees of pointing accuracy). One possibility for 
this capability is a low-power, non-complex safe mode for a more complex ACS. This 
research is helpful though to find the limits of technology necessary to make this a viable 
alternative method. The possibility of a low cost solution to attitude control that could 
salvage a mission after complex mechanical failures occur would be a substantial 


motivation for research to find such a method. 


The NPSATI satellite will use an active magnetic torque rod system with a three- 
axis magnetometer for attitude determination to maintain 3-axis stabilization with a 
slightly gravity gradient friendly structure. This thesis will examine the performance of 
three combinations of programs and simulation models that may help start providing the 
evidence that such a system is feasible and viable for the NPSAT1 satellite. Graphs and 


tables will be included to augment the presentation of the analysis data in this thesis. 
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Hi. BACKGROUND AND HISTORICAL PERSPECTIVE 


NPSatl1 is the second satellite in the Naval Postgraduate School’s Small Satellite 
Design Program. The purpose of the Small Satellite Design Program is twofold. First, it 
provides hands-on education for officer students in the disciplines of satellite design, 
system integration, space systems operations, and software development. Second, it 


provides students with dedicated, on-orbit assets for experimentation and testing. 


NPSatl is an incremental improvement over the Postgraduate School’s first small 
satellite, PANSAT. NPSatl is a three-axis stabilized spacecraft equipped primarily with 
Commercial-Off-The-Shelf (COTS) hardware. Because of its three-axis stabilization, it 
is anticipated that the spacecraft could support a basic imaging capability with a COTS 


digital camera. 


Several spacecraft components have already been acquired or built, and several 
experiments or payloads have already been envisioned for the mission. The solar cells 
and basic structure of the spacecraft were inherited from a cancelled Navy program at no 
cost to the Naval Postgraduate School. The primary spacecraft controller, a modified 
COTS 80386 PC-104 board, was already designed and built. The Digital Camera and its 
interface board had been purchased and tested with the PC-104 Central Processor. The 
Naval Research Laboratory has committed to flying their Coherent Electromagnetic 
Radio Tomography (CERTO) experiment onboard. Dr. Alan Ross, Navy TENCAP Chair 
in the Space Systems Academic Group, has indicated a desire to test his Triple 


Redundant Modular Computer experiment on board NPSAT1 as well. 


The basic orbit profile for NPSAT1 was to be a circular Low Earth Orbit (LEO) 
between 400 and 800 kilometers, inclined at between 40 and 80 degrees with a mission 


life of 18 to 24 months. 


The NPSat1 Attitude Determination and Control System (ADCS) was designed to 
provide three-axis stabilization. The spacecraft is gravity gradient friendly with the 
possibility of an optional telescoping boom for added stability. Active attitude control is 


achieved using three magnetic torque rods. Attitude determination software receives 


inputs from a three-axis magnetometer to determine spacecraft attitude and rates. Figure 
2.1 shows multiple views of the NPSAT1 spacecraft and the planned solar array layout. 
Three-axis stabilization provides a stable spacecraft to support the imagery and CERTO 
Missions on NPSatl. The cone is at the top of the spacecraft and the nadir pointing face 


is the bottom. 





Figure 2.1. NPSATI 


The NPSAT1 spacecraft will be deployed from an Air Force Delta class 
expendable launch vehicle. The spacecraft will be deployed from an EELV secondary 
payload adapter (ESPA) ring. Figure 2.2 shows the NPSAT1 spacecraft as it will be 
configured on the ESPA ring. 





Figure 2.2. NPSAT1 attached to EELV payload adapter (ESPA) 


THIS PAGE INTENTIONALLY LEFT BLANK 


Hl. STATEMENT OF PURPOSE 


A. OBJECTIVES OF RESEARCH 


Initially, two objectives were sought, including sensitivity analysis of a magnetic 
control law program, and after a suggestion from the Naval Research Laboratory to try a 
magnetic rate (Bdot) control law [Ref. 1], a benefit analysis of including a magnetic rate 
control law for initial angular rate reduction of the NPSAT1 spacecraft. In the course of 
conducting the analysis, two secondary objectives were identified. Analysis was then to 
include the difference in system performance between switching from Bdot control to 
nadir-pointing control after a time interval or after acquisition of a specified reduced 
angular rate. The remaining secondary objective was to explore the benefits of mode 
switching before or after priming the systems reduced order estimator with historical 


data. 


The analysis was conducted to determine the best configuration of the existing 
proposed attitude control systems for use on the NPSAT1 spacecraft. Two similar 
attitude control systems were examined in three different configurations. The first 
system consisted of two MATLAB programs and one SIMULINK model without a Bdot 
control loop. The second system consisted of one MATLAB program and one 
SIMULINK model with a timed hand-off from an integrated Bdot control loop. 
Additionally, the second system’s reduced estimator received priming data prior to hand- 
off. The third configuration consisted of an independent Bdot control program which 
executed a control hand-off to the first SIMULINK model without priming the reduced 
estimator, after a predetermined, reduced, combined angular rate was achieved. Table 


3.A.1 shows the programs included in each system configuration. 


The research focused on three objectives. These included: initial condition 
sensitivity to all of the three program configurations, the difference between including the 
Bdot or no Bdot, and the difference between a timed Bdot hand-off or an angular rate 


based hand-off. 


SIMULINK model MATLAB system programs 
LAB 17.mdl Lab15MOIData.m 
P= Lab 6 Data.m 
NPSAT1ACSINTRYW.mdl |NPSAT1ACSData.m 
Lab17.mdl Lab15MO0IData.m 


"Lab17 hand-off model" Labi6Data.m 
nps.m 


Table 3.A.1. Configuration models and programs 








B. APPROACH 


Each of the two SIMULINK models was designed to generate control torque 
commands to three magnetic torque rods for active attitude control. Using initial 
spacecraft ephemeris data and a spherical harmonic model (8"" order), the systems 
estimate what the Earth’s magnetic field should be at its current location. The system 
then computes the cross product of the expected magnetic field with the actual local 
magnetic field measured by a magnetometer. A cross product of the two magnetic 
vectors is used as an approximation of a magnetic control law to produce the first three 
system states (approximate euler angles). A reduced estimator is then used to 
approximate the other three system states (angular rates). These six states are fed to an 


actuator control law to generate torques. 


The baseline model for the comparison was the SIMULINK model 
‘NPSATACSIntRvw.mdl’ and its associated m-file, ‘NPSAT1ACSDATA.m’ which 
performed a timed Bdot control hand-off as previously described. [Ref. 2] This model is 
hereafter referred to as the “combo” model since the Bdot control loop is integrated 
within the SIMULINK model. The combo model, whose block diagram is shown in 
Figure 3.B.1, was used for the analysis baseline after analysis indicated that configuration 
exhibited the lowest power consumption of the three configurations. The three 


configurations of systems were evaluated for time response and power consumption. 


B=Measured magnetic field vector 
b=Calculated magnetic field vector 





Figure 3.B.1. Block Diagram for “Combo” model 


The programs and models for each of the three configurations are included in 


Appendix A. 
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IV. PROGRAM SETUP AND PROCEDURE 


A. INITIAL CONDITIONS AND ASSUMPTIONS 


In order for the attitude control system for the NPSAT1 spacecraft to be viable, it 
must be able to handle the conditions immediately after deployment from the launch 
vehicle, and handle the maintenance of the spacecraft’s operational pointing 
requirements. Initial conditions and assumptions were established. A representative 
number of test cases were run to establish the capability of the system to respond to those 
conditions. Each configuration was compared against the others to establish the best 
candidate for the final attitude control system. The primary measures of effectiveness 


were time to reach steady state and system power consumption. 


The assumptions and initial conditions used during the simulations were chosen to 
approximate as closely as possible, the expected final characteristics of the satellite and 
its anticipated orbit. The minimum allowable altitude due to aerodynamic drag and 
aerodynamic disturbance torques for the spacecraft is 400km. The expected orbit is 
circular at 600km. As Figure 4.A.1 shows, the three 30 Am“2 torque rods can easily 
counter solar, gravity gradient, or aerodynamic torques encountered at 600km. The 
generated torques are greater than an order of magnitude more powerful than the 


anticipated disturbance forces. 
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Torque vs Altitude 
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Figure 4.A.1. Torque vs. Altitude. 


Initial satellite angular rates were based on the expected possible tip-off rates 
encountered by the launch vehicle at the moment of deployment. It was assumed that +/- 
0.1 rad/sec was a conservative range of possible rates for spacecraft pitch and roll. Since 
the yaw rate induced by the launcher to the spacecraft is not expected to be as great as the 
pitch and roll rates, it was given a range of +/- 0.01 rad/sec. For simplification and 
standardization of initial conditions between tests, the initial pitch and roll rates were 
assigned either a “high rate” magnitude (+ 0.1 rad/sec or —0.1 rad/sec) or a “low rate” 
magnitude (+0.01 rad/sec or —0.01 rad/sec). The initial yaw rate was assigned a positive 
or negative “low rate” only. Considering all the possible combinations of the “high” and 
“low” angular rates resulted in 32 different, testable, initial angular rate combinations. 
For simplification of notation, the angular rates were assigned a three letter combination 
of either capital or lower case N’s or P’s, representing the magnitude of each rate as well 
as its sign (negative or positive). Example: nNp represents an initial pitch, roll and yaw 


rate of —0.01 rad/sec, -0.1 rad/sec and +0.01 rad/sec respectively. 


Power consumption, as well as time response analysis of the various models to 


each of the selected test cases required the capture of data from 11 diferent program 
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variables. The variables included the output torques for each of the three torque rods, a 
sum of the absolute values of the three rods’ output, the Euler angles , 0, and y, the 


angular rates ddot, Odot, and wdot, and time. 


B. EXPERIMENTAL PROCEDURE 


1. Correlation Between Initial Conditions and Test Cases 


The Bdot control law program ‘nps.m’ and the combo SIMULINK model were 
each evaluated for all 32 expected, representative, initial angular rate combinations. The 
programs’ time responses represented the time each program required to reduce the 
angular rates to a combined rate of less than 1.5 degrees per second. The resulting time 
responses seemingly had no correlation with any particular aspect of the initial 
conditions. Table 4.B.1.1 shows the lack of apparent correlation between the time 
responses for the Bdot control law program ‘nps.m’ and the combo model and the initial 


conditions. 


Time req'd : 
Bdot [sec 





——7 a 
int Fates |i int Fates int Fates | Bdot (sec 
INPn | Inen [6900] [Nop | 60 
npn | nen [e900] [Pan | 14000 
Np | Iwnip [11000] [nm [22 
Pan | Pan [1100] [pp | __¢700 
Inn | Inn [3660] [PPp | 11000 
win | Inn [11000] [nnn | 3000 
lpm | lpm [380] [PPn [11200 





Table 4.B.1.1. Time to Achieve System Steady State. 


Except where all of the initial rate magnitudes were “low” (0.01 rad/sec), in 
which case response times tended to also be relatively low, the data did not support 


linking higher initial angular rates with longer system settling times, as was expected. 
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Neither did the results suggest any correlation between sign order of the initial angular 
rates with a higher or lower system settling time. Instead, the resulting data seemed to 
suggest that the ‘nps.m’ program had certain initial positive or negative angular rate 
combinations that the program responded to poorly with respect to system settling times. 
Many combinations whose sign order matched, yet magnitudes varied, produced similar 
system steady state time responses, regardless of the difference in individual angular rate 


magnitudes. 


The data did not favor a spread of test cases based on initial angular rate 
combinations as a means to guarantee a representative spread of system performance. 
Instead, using the data collected by running the combo model simulation (configuration 
2), the 32 combinations of initial angular rates were sorted in order of their system steady 
state settling times. The combo model was chosen as the standard of performance to 
measure all other programs or combinations of programs against to determine the optimal 
system configuration for the full range of expected initial conditions. Two pieces of data 
of specific interest in the analysis was the total system torque rod output, used to 
calculate an energy index, and the time to reach the system steady state for each test case. 
Figure 4.B.1.1 shows an early plot of instantaneous and cumulative torque rod power. 


Note that the torque rods chosen for the NPSAT1 spacecraft are limited to +/- 30 Am‘’2 
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Total Power = 0.53037 VV-hrs 
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Figure 4.B.1.1. Three Torque Rod Output vs. Time. 


Nine test cases were chosen to represent a spread of settling times demonstrated 


by the combo model, from the lowest demonstrated settling time to the highest. The nine 


99 66. 99 66 


sets of initial angular rates included “ppn,” ““NPn,” “NNp,” “pPn,” “Pnp,” “nnn,” “pnp,” 
“nNp,” and “Ppp.” 
2. Run Program, Capture Data 


There were three configurations of programs and models that were run for each of 
the nine test cases. For each iteration, the data was recorded and analyzed to determine 
the power required by the system to reach the desired system steady state. The power 
required is represented in the form of an energy index. The torque rod output is 
measured in Am‘2. If the output is integrated over the length of time to reach system 
steady state, the result is an “energy index” (Am%2 * sec) that was used for comparison 
against the other three programs and methods that were run with the same initial 


conditions. 


In order to capture the data variables of interest from the simulations and program 
runs, the programs and simulation models were modified to include output blocks to 
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generate MATLAB ‘.mat’ files which included the desired data, or the programs were 
modified to generate new variables which could be saved from the MATLAB work 
space. The programs ‘nps.m’ and the combo model were each run for all 32 


combinations of initial angular rates. 


The Bdot control law program ‘nps.m’ required modifications of the MATLAB 
source code to record the 11 variables from each program run. Additional data 
manipulation and plotting programs were created for further data analysis and to plot the 
analysis results. The program ‘grabdataset.m’ was written for execution following every 
different iteration of the ‘nps.m’ program to clear the MATLAB workspace after copying 
all desired variables into an 11 element matrix. Each of the variables was sampled at a 
pre-assigned change in time or dt (sec) until the simulation reached its chosen final 
combined rates goal. The value chosen for all of the simulations was dt = 2.0 seconds. 
The plotting programs ‘plotresult2.m’ and ‘rodtotal_vs_time.m’ were used to provide 
graphical representation of the data. The source code for the programs ‘grabdataset.m,’ 


‘plotresult2.m’ and ‘rodtotal_vs_time.m’ are provided in Appendix A. 


Both SIMULINK models were modified to capture the same 11 variables that 
were recorded from the ‘nps.m’ program runs. Each SIMULINK model was run with the 
same initial conditions used in the iterations that were run with the Bdot control law 


program. The output blocks added to the simulation models are shown in Figure 4.B.2.1. 


3 rod torques Euler angles = Angular rates 


















410utput_Vars.mat 





To File 





thd 
psd 
phi 
theta 
psi 








Figure 4.B.2.1. Simulink Output Data Capture Blocks. 
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The output ‘.mat’ files were imported into the MATLAB work space for data 
analysis and plotting using the programs ‘grabdataset.m,’ ‘create_variables.m,’ 


‘rodtotal_vs_time.m’ and ‘plotresult2.m,’ which are provided in Appendix A. 


The program ‘grabdataset.m’ was written for manipulating data captured from a 
SIMULINK simulation run. The program measured and plotted the instantaneous and 
cumulative total of commanded rod torques and the time required to achieve a 3-axis 
steady state pointing accuracy of less than 10 degrees. This program also solved and 
plotted the area under the instantaneous rod torque total curve, producing an energy index 


for comparison with other runs using different initial conditions. 


Eleven variables were recorded in a matrix as a ‘.mat file’ from a simulation using 
one of two modified SIMULINK simulations. The program picked off the following 
desired variables: the torque rod outputs in Am“2 as ‘rod1’, ‘rod2’, and ‘rod3.’ The total 
commanded torque from the three rods combined was ‘rodtot.’ The angular rates were 
recorded as ‘Phd,’ ‘Thd’ and ‘Psd,’ and the Euler angles were recorded as ‘Ph,’ ‘Th’ and 
‘Ps.’ 


The Bdot MATLAB program output data and SIMULINK system models’ output 
data needed to be merged into a common format for data analysis. The program ‘create 
variables.m’ was written for capturing data from the NRL program ‘nps.m.’ The 
program captured 11 variables for analysis and plotting purposes. Many of the data 
parameters were manually entered into the program code from the MATLAB workspace. 
Future analysis will benefit from the development of an automated data capture program 
that combines all of the steps used to manipulate the data for this analysis. The resulting 
data was saved from the MATLAB workspace to a data folder. The ‘.mat’ file’s name 


reflected the initial conditions used to produce the data. 


EXAMPLE FILE NAME: ‘bdot_pNp_casel.mat.’ In this example, pNp stands 
for ‘Phd’ =.01 rad/s, ‘Thd’ =-.1 rad/s, ‘Psd’ =.01 rad/s. 


The program ‘Rodtotal_vs_Time.m’ was written to plot simulated torque rod data 
for analysis of simulations supporting NPSSAT1. This program captured a row of data 


from a saved 11 row matrix. The matrix was produced from one of two modified 
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SIMULINK simulations. The program grabbed torque rod totals and plotted them vs. 
time. 


3. Run Program Handoff Method, Capture Data 


The program ‘NPSSATITESTCASE.m’ was written for manipulating data 
captured from MATLAB programs and SIMULINK runs that were run in the hand-off 


mode. Since each program ran independently, their data had to be combined for analysis. 


This program computed and plotted the results of the two different simulations’ 
data. The Bdot control law program ‘nps.m’ was used to arrest angular rotation rates to 
the predetermined rate goal (less than or equal to 1.5 degrees/sec). Then, the final states 
of that simulation were used for the initial conditions of the second program, which 
further arrested the rotation rates and achieved a predetermined three-axis nadir pointing 
accuracy. Data from the program ‘nps.m,’ was used to provide the initial conditions for 
each of two modified SIMULINK simulations. The data stored as a result of running one 
of the two follow-on simulations were loaded into the MATLAB work space, with the 
stored data from the Bdot program run. The ‘NPSSAT1TESTCASE.m’ program 
essentially combined two other programs, ‘grabdataset.m’ and ‘create_variables.m.’ 

After loading the desired ‘.mat’ files containing data matrices into the MATLAB 
workspace, the matrix sizes were used to edit the program’s array sizes. The data are 
saved as a single ‘.mat’ file. 

The ‘NPSSATITESTCASE.m’ program captured torque rod total and time 
variables for analysis and plotting purposes. The program then measured and plotted the 
instantaneous and cumulative total of commanded rod torques and the time required to 
achieve a 3-axis steady state pointing accuracy of less than 10 degrees. This program 
computed and plotted the area under the instantaneous rod total curve to produce an 
energy index for comparison with other runs using different initial conditions. 


4. Run Plotting Programs 


Once the four programs or combinations of programs were run for each of the 
nine test cases, the ability to use and manipulate the data with calculations was relatively 
straight forward. In order to manipulate the data from both the independent Bdot 
program and its second handoff program, the data had to be merged together at the point 
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of handoff, which occurred at the instant the Bdot program reached its assigned combined 
angular rate threshold. The resultant data was then used as if it were merely the results of 


one of the individual programs. 
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Figure 4.B.4.1. Angular Rates vs. Time (pPn). 


An example of the plotting program “plotresult2.m’ is shown in Figure 4.B.4.1. 
The initial angular rates were ddot = +0.01 rad/sec, @dot = +0.1 rad/sec, and wdot = -0.01 
rad/sec (pPn). The program ‘nps.m’ was assigned a combined angular rate goal of 1.5 
deg/sec. The program arrested the angular rates until the sum of the absolute values of 


the three angular rates dropped below the assigned threshold. 


Figure 4.B.4.2 shows the results of configuration 3 (on the left side of the figure), 
the Bdot control handoff method to the ‘Lab17.mdl’ for the test cases “pPn” (dot = 
+0.01 rad/sec, Odot = +0.1 rad/sec, and wdot = -0.01 rad/sec) and “NNp.” (odot = -0.1 
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rad/sec, O@dot = -0.1 rad/sec, and wdot = +0.01 rad/sec). In this method, the SIMULINK 
model ‘Lab17’ has not been activated until the Bdot program ‘nps.m’ has finished 


arresting the angular rates. 
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Figure 4.B.4.2. Bdot ‘nps.m’ handoff to ‘Lab17.mdl’. 


The point at which the handoff occurred is quite distinguishable in these cases, 
occurring at 5,000 seconds and 1,500 seconds respectively. The graphical representation 
of the data in this case was instrumental in recognizing a difference in the smoothness of 
transitions between different handoff models. Using the same initial conditions for 
configuration 2, the internal system Bdot control handoff to the estimator in the combo 
SIMULINK model shows a much smoother transition at the handoff point. As Figure 
4.B.4.2 indicates, the sudden spike in required energy from the torque rods does not 


occur. This observation proved to be true for many of the handoff cases. 


20 


While attempting to explain why there was a difference between the 
configurations, several insights were gained into the operations of the models. The 
combo model, essentially, is a Bdot loop operating for a period of time within the model, 
then handing off the control to the reduced order estimator. The Kalman filter within the 
seemingly idle reduced order estimator loop of the model was in fact gathering magnetic 
field crossed with magnetometer information, while the Bdot loop was arresting the 
initial angular rates. Therefore, when the control was handed off by the Bdot loop to the 
reduced order estimator, the filter already had better knowledge than if the system started 
gathering magnetic field data at the handoff point. The right side of Figure 4.B.4.2 shows 
the handoff which includes the “better knowledge.” 
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V. RESULTS 


A. OVERVIEW 


The progression of analysis objectives began with sensitivity analysis of the 
existing ACS model, which had no Bdot control law. While that analysis was being 
conducted, Naval Research Laboratory personnel recommended using a magnetic rate 
(Bdot) control law to reduce the satellite angular rates to a lower “manageable” level 
prior to using the reduced estimator models. [Ref. 1] The objective of showing the 
difference in performance between systems with or without the Bdot control law was 
started. It was expected that the systems would acquire steady state nadir pointing faster 
if the initial angular rates were reduced faster. The Bdot program ‘nps.m’ was analyzed 
for 32 different sets of initial conditions and reduced the rates much faster than the 
reduced estimator control system without Bdot. The performance measurements were 
based on the time required for the system to achieve nadir poining accuracies of +/- 10 


degrees. 


After several evolutionary steps of the ACS system, the combo model was 
produced, which incorporated the Bdot control law within the system model. The power 
consumption results out-performed all previous system versions. The decision to use the 
combo model as a benchmark for performance was made. The combo model, 
configuration 2 was analyzed for all 32 initial test case combinations for time response. 
Curiously, the time responses for the combo model compared to the other configurations 
did not have the best performance. It was discovered that a difference in time response 
was achieved if the Bdot hand-off was conducted after a predetermined, reduced angular 
rate was achieved instead of after a predetermined time interval. The decision to conduct 
analysis in support of three objectives was then made. The three objectives included: 
initial condition sensitivity to all of three program configurations, the difference between 
Bdot and no Bdot, and the difference between timed Bdot hand-off and angular rate 
based Bdot hand-off. 


23 


B. BDOT RESULTS 


The results of the Bdot program runs offered no insight into the correlation 
between initial conditions and system time response. Sorting the data by time required to 
reach the combined angular rate goal of <=1.5 degrees/sec, had no correlation with a 
similar sorting of any of the other systems’ results. There is no doubt that the Bdot 
control law program, “nps.m,’ by itself, arrested the angular rates faster than any other 


program or model. (See Figure 5.B.1) 
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Figure 5.B.1. Bdot Time to Reduce Rates to <= 1.5 deg/sec Combined. 


Unfortunately the program does not offer any three-axis pointing capability. 
There was, however, substantial time-savings, often by a factor of ten or more over the 
ability of the other programs to achieve partial acquisition by angular rate reduction. 
This performance led to the decision to analyze the affect of using the Bdot program 
independently until a specified combined angular rate goal was achieved. The final 
conditions, including angular rates and Euler position angles were handed off to the 


original ‘Lab17.mdl’ as initial conditions. (configuration 3). 


The Bdot control law more efficiently arrested the angular rates, compared to 


configuration 1, the system without the magnetic control law (‘Lab17.mdl’). The Bdot 
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program also appeared to help the other systems most efficiently when the Bdot program 
was run to a specific rate goal, instead of for a specified time. Once the control was 
handed off to the other programs or models, the energy levels and commanded torques 


required to resume the control effort were minimal. 


C. CONFIGURATION 1 RESULTS 


After conducting the nine test cases, configuration 1 seemingly had the poorest 
performance of the three configurations. Compared to the other two configurations, to 
achieve pointing accuracies of +/- 10 degrees, configuration 1 appeared to only perform 
well in the middle ground between the extremes of initial conditions (using the combo 
model as a bench mark). The best cases and worst cases for configuration 2 (the 
extremes), however, nearly always corresponded to poor performance from this 


configuration. 


The energy indexes for this configuration tended to exceed the worst performance 
of the other two methods. The gains selected for this model, however, have been 
designed to be adequate for the entire range of initial conditions. It is quite likely that if a 
Bdot control law was added to this model, and the Bdot program maintained its control 
until a minimal residual combined angular rate existed, then the gains for the 
‘Lab17.mdl’ could be optimized to accommodate the low initial rate conditions. Table 
5.C.1 shows the nine case studies conducted. The order of the initial conditions is in 
increasing system time order for the reference combo model. Figures 5.C.1 through 
5.C.18 show the instantaneous and cumulative power consumption results for 


configuration 1. 
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Table 5.C.1. Configuration | (‘Lab17.md1’) test results. 
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Figure 5.C.1. “ppn” Configuration 1: Cumulative Rod Total vs. Time 
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“ppn” Configuration 1: Instantaneous Rod Total vs. Time 


Figure 5.C.2. 
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Figure 5.C.3. “NPn” Configuration 1: Cumulative Rod Total vs. Time 
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Figure 5.C.4. “NPn” Configuration 1: Instantaneous Rod Total vs. Time 


28 





1 
1 
1 
' 
1 
' 
L 
1 
' 
1 
1 
1 
' 
' 
a 
1 
1 
1 
1 
1 
1 
' 
r 
1 
1 
1 
1 
1 
' 
' 
T 
' 
1 
1 
1 
' 
' 
L 
1 
1 
1 
1 
1 
1 
' 
7 
1 
1 
1 
1 
1 
' 
' 
t 
' 
1 
1 
1 
1 
' 
' 
t 
1 
1 
' 
1 
' 
L 
1 
1 
1 
1 
1 
' 


eS ee ee ee ee ee eee ee ee ee ee 


re es te ee ee et ee ee ee ee ee eee eee 





Figure 5.C.5. “NNp” Configuration 1: Cumulative Rod Total vs. Time 
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Figure 5.C.6. “NNp” Configuration 1: Instantaneous Rod Total vs. Time 
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Figure 5.C.7. “pPn” Configuration 1: Cumulative Rod Total vs. Time 
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Figure 5.C.8. “pPn” Configuration 1: Instantaneous Rod Total vs. Time 
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Figure 5.C.9. “Pnp” Configuration 1: Cumulative Rod Total vs. Time 
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Figure 5.C.10. “Pnp” Configuration 1: Instantaneous Rod Total vs. Time 
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Figure 5.C.11. “nnn” Configuration 1: Cumulative Rod Total vs. Time 





Figure 5.C.12. “nnn” Configuration 1: Instantaneous Rod Total vs. Time 
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pnp(Lab17} area Rod Total 98606.254 in 170000 sec 
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Figure 5.C.13. “pnp” Configuration 1: Cumulative Rod Total vs. Time 
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pnp” Configuration 1: Instantaneous Rod Total vs. Time 
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Figure 5.C.14. 


33 


1 
1 
1 
1 
1 
1 
d 
1 
1 
1 
1 
' 
1 
1 
4 
1 
1 
1 
1 
1 
1 
1 
4 
1 
1 
1 
1 
1 
1 
1 
4 
1 
1 
1 
1 
1 
1 
1 





Figure 5.C.15. “nNp” Configuration 1: Cumulative Rod Total vs. Time 





“nNp” Configuration 1: Instantaneous Rod Total vs. Time 


Figure 5.C.16. 
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Figure 5.C.17. “Ppp” Configuration 1: Cumulative Rod Total vs. Time 
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Figure 5.C.18. “Ppp” Configuration 1: Instantaneous Rod Total vs. Time 
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D. CONFIGURATION 2 RESULTS 


The combo model’s results became the baseline reference for all data comparison. 
The power required to achieve steady state was consistently lower than any other model 
or program, regardless of initial satellite angular rates. The advantage of lower required 
power levels, however, was somewhat diminished by the apparent trade-off with 
increased settling times. The majority of the other models and programs consistently 
performed better than the configuration 2 model with respect to time to achieve steady 
state. Nevertheless, time to achieve steady state has been considered of less overall 


importance than power use for the initial acquisition of orbit attitude. 


All 32 initial angular rate test cases were analyzed using the combo model. There 
was, like the results of the Bdot program runs, no apparent correlation between sign or 
magnitude combinations of the initial angular rates and the resulting time or power 
consumption required to achieve a steady state condition. There appeared, however, to 
be some correlation between power consumption and time to achieve steady state. The 
initial power use for “high” satellite angular rates and “low” satellite angular rates is 
comparable between test cases, regardless of which angular rate had the highest 
magnitude. All configuration 2 test cases had comparable power consumption versus 
time, once the initial higher rates were “knocked down” by the system. Therefore, the 
test cases that took the longest time to achieve a steady state condition, correspondingly 
also had the highest power consumption levels. Figures 5.D.1 through 5.D.18 show the 


instantaneous and cumulative power consumption results for configuration 2. 
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Figure 5.D.1. “ppn” Configuration 2: Cumulative Rod Total vs. Time 





Figure 5.D.2. “ppn” Configuration 2: Instantaneous Rod Total vs. Time 
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Figure 5.D.3. “NPn” Configuration 2: Cumulative Rod Total vs. Time 





Figure 5.D.4. “NPn” Configuration 2: Instantaneous Rod Total vs. Time 


38 


1 
1 
1 
1 
' 
' 
1 
1 
1 
+ 
1 
' 
' 
1 
' 
1 
1 
1 
1 
7 
' 
' 
1 
1 
1 
1 
1 
' 
4 
1 
1 
1 
1 
1 
' 
' 
' 
' 
4 
1 
1 
1 
' 
' 
' 
1 
1 
1 
7 
1 
' 
' 
1 
' 
1 
1 
1 
1 
i 
' 
' 
1 
1 
1 





Figure 5.D.5. “NNp” Configuration 2: Cumulative Rod Total vs. Time 
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Figure 5.D.6. “NNp” Configuration 2: Instantaneous Rod Total vs. Time 
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Figure 5.D.7. “pPn” Configuration 2: Cumulative Rod Total vs. Time 
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Figure 5.D.8. “pPn” Configuration 2: Instantaneous Rod Total vs. Time 


40 


Pnp(Combo) area Rod Total 77881.8498 in 100000 sec 


eee ee eee eee 


ae 
7 --- 

2 

1 






Time (sec) 


Figure 5.D.9. “Pnp” Configuration 2: Cumulative Rod Total vs. Time 
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Figure 5.D.10. “Pnp 
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Figure 5.D.11. “nnn” Configuration 2: Cumulative Rod Total vs. Time 





nnn” Configuration 2: Instantaneous Rod Total vs. Time 


ce 


Figure 5.D.12. 
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Figure 5.D.13. “pnp” Configuration 2: Cumulative Rod Total vs. Time 
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Figure 5.D.14. “pnp” Configuration 2: Instantaneous Rod Total vs. Time 
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Figure 5.D.15. “nNp” Configuration 2: Cumulative Rod Total vs. Time 
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nNp” Configuration 2: Instantaneous Rod Total vs. Time 
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Figure 5.D.16. 
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Figure 5.D.17. “Ppp” Configuration 2: Cumulative Rod Total vs. Time 





Figure 5.D.18. “Ppp” Configuration 2: Instantaneous Rod Total vs. Time 
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E. CONFIGURATION 3 RESULTS 


When the Bdot program was run independently, then handed-off to the LAB17 
model, using the final conditions from the Bdot test case as the initial conditions for the 
LAB17 test case, the following results occurred: Compared to the nine combo model test 
cases chosen to represent the full spread of outcome possibilities, the LAB17 hand-off 
method test cases tended not to perform as well compared to the fastest half of the combo 
method test cases. It performed better, though, than the slowest half of the combo 
method test cases. Performance levels, mentioned above, include both power 
consumption and time response to test case initial conditions. The results also further 
amplified the finding that there was no apparent correlation between sign or magnitude of 
the test case initial angular rates and the performance of the system in power 
consumption or time response. Figures 5.E.1 through 5.E.2 show the instantaneous and 


cumulative power consumption results for configuration 3. 


Cases “NPn,” “NNp,” “pPn,” “nnn,” and “Ppp” clearly illustrate the power spike 
occurring at the hand-off point between the Bdot and Lab17 models. This is most likely 


due to the poor initial estimates within the Kalman Filter. 
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ppnhandoff(Lab1?) area Rod Total 86653.1201 in 90016 sec 
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Figure 5.E.1. “ppn” Configuration 3: Cumulative Rod Total vs. Time 
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Figure 5.E.2. “ppn” Configuration 3: Instantaneous Rod Total vs. Time 
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Figure 5.E.3. “NPn” Configuration 3: Cumulative Rod Total vs. Time 
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Figure 5.E.4. “NPn” Configuration 3: Instantaneous Rod Total vs. Time 
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Figure 5.E.5. “NNp” Configuration 3: Cumulative Rod Total vs. Time 





Figure 5.E.6. “NNp” Configuration 3: Instantaneous Rod Total vs. Time 
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Figure 5.E.7. “pPn” Configuration 3: Cumulative Rod Total vs. Time 





Figure 5.E.8. “pPn” Configuration 3: Instantaneous Rod Total vs. Time 
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Figure 5.E.9. “Pnp” Configuration 3: Cumulative Rod Total vs. Time 
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Configuration 3: Instantaneous Rod Total vs. Time 
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Figure 5.E.10. 
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Figure 5.E.11. “nnn” Configuration 3: Cumulative Rod Total vs. Time 
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Figure 5.E.12. “nnn” Configuration 3: Instantaneous Rod Total vs. Time 


52 


eaeecaseseciviosea sce saae 





Figure 5.E.13. “pnp” Configuration 3: Cumulative Rod Total vs. Time 
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Configuration 3: Instantaneous Rod Total vs. Time 
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Figure 5.E.14. 
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Figure 5.E.15. “nNp” Configuration 3: Cumulative Rod Total vs. Time 
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Figure 5.E.16. “nNp” Configuration 3: Instantaneous Rod Total vs. Time 
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Figure 5.E.17. “Ppp” Configuration 3: Cumulative Rod Total vs. Time 
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Figure 5.E.18. “Ppp” Configuration 3: Instantaneous Rod Total vs. Time 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 


Table 6.A.1 summarizes the comparative results of each of the three attitude 
control methods analyzed in this study. The reference combo method was chosen 
primarily for its relatively superior performance in power consumption required to reach 
a steady state. The nine comparison test cases were chosen to represent the entire 


performance range of the combo method test cases. 


The most obvious conclusions that can be obtained from the analysis of these 
control methods and resulting data are as follows: The combo model, configuration 2 
generally consumed the lowest power levels of any other method, but required the 
longest settling times for medium to high initial angular rates. The control methods that 
included the Bdot control law (configuration 2 and 3) quickly reduced the satellite 
angular rates to a level that thence required only minimal power consumption to achieve 
steady state pointing accuracies of less than ten degrees. Time to achieve system steady 
state was achieved more quickly and with less power consumption for cases with higher 
initial angular rates when the hand-off from a Bdot control law occurred after a threshold 
combined angular rate was achieved instead of performing the hand-off after a 
predetermined time interval. The configuration 2 model used a timed hand-off from its 
internal Bdot law control loop, and the configuration 3 model handed-off from a separate 
Bdot program to the various other models after a threshold combined angular rate was 


achieved. 


The first of the three main analysis objectives was to examine sensitivity of the 
three configurations to initial conditions. All three configurations performed well for all 
the test case initial conditions and were viable options for an ACS. The second objective 
was to compare the differences between using Bdot or no Bdot within a configuration. 
Although there was questionable improvement in system time responses between the 
configurations, the power consumption results indicated that both configurations 2 and 3, 
which used a Bdot control loop, required less power to achieve steady state than 


configuration 1, which did not include Bdot. The third objective was to explore the 


af 


difference between using a timed Bdot hand-off or an angular rate based Bdot hand-off. 
Configuration 2, which used a timed hand-off, had better power consumption 
performance, but better time response results than configuration 3, which used the 
angular rate based hand-off. 
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Table 6.A.1. Test Case Table Data. 
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Figure 6.4.2 shows common system power consumption response characteristics 
that can be gleaned from the visual comparison of the three configurations. Magnetic 
torque rod power consumption was lower, and high angular rates were generally arrested 
faster, when a Bdot control law was used to reduce initial angular rates to a more- 


manageable “low” level. When the Bdot control law was used to achieve a sufficiently 
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low threshold combined angular rate, as opposed to a timed hand-off from a Bdot control 


law, the power spike that often occurs at the hand-off point was eliminated. 
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Figure 6.A.2. “pPn” Test Case. 








B. RECOMMENDATIONS 


The procedure for conducting an analysis of the handoff from ‘nps.m,’ which was 
the Bdot control law program to another program was cumbersome and time-consuming. 
If further analysis is deemed appropriate to determine whether a Bdot control law should 
be used for NPSAT1, a combined MATLAB program incorporating each of the smaller 


steps should be developed to more fully automate the data acquisition and analysis. 


It appears, howerer, that there is clear evidence that the Bdot control law is a 
powerful and efficient method to arrest initial angular rates, regardless of their 
magnitude. Should it be decided to go forward with some implementation of a Bdot 


control law, it is recommended that the Bdot control law be a selectable control loop 
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within a combined program or system model. It is recommended that the hand-off that 
will occur from the Bdot control law loop should be done after a threshold combined 
angular rate is achieved, and not after a predetermined time interval. The time for the 
handoff models to begin steady state pointing accuracy acquisition after the initial 
arresting of the launch vehicle tip-off rates or a loss of attitude control can be minimized 
by handing off control as soon as the benefit of the Bdot control law ceases. The gains 
within the control loop that will receive a hand-off from the Bdot control law loop can be 
optimized for very low angular rates, since the Bdot control law would handle any 


initially higher angular rates. 


All of the control methods analyzed appeared to be valid control methods to 
achieve three-axis attitude stabilization using only magnetic torquers and magnetometers 
for active control. The most efficient control method analyzed appeared to be 
configuration 2. This method was achieved by simply modifying the configuration 1 
model to perform the Bdot hand-off based on achievement of a predetermined combined 
angular rate as opposed to a predetermined run time. This method also appears to 
eliminate the power spike that seemingly occurs when the angular rates have not been 


sufficiently arrested prior to hand-off. 
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APPENDIX A. PROGRAMS 


NRL BDOT CONTROL LAW (NPS.M) 


% Matlab simulation files courtesy of Dr. Glenn Creamer, NRL Code 8231 
% Main program (nps.m) 


clear time wl w2 w3 m1 m2 m3 b1 b2 b3 

clear ugg] ugg2 ugg3 umagl umag2 umag3 

clear yaw_angle pitch_angle roll_angle libration 

global earth_radius satellite inertia inv_satellite_inertia dt 


% output storage parameters 
count=0; 

k=0; 

store_count=3; 


% conversion factors 
d2r=pi/180; 
sqrt2=sqrt(2); 


% Earth parameters 
earth_radius=6378.1363; 
earth_rate=7.272205217e-05; 
earth_gravity_constant=3.986004415e+05; 


% orbit parameters 


inc=35*d2r; YW*eREK LY PDATE DATA ITEM 
ra=0*d2r; 
altitude=600; YW*eREX UPDATE DATA ITEM 


radius=earth_radius+altitude; 


% orbit rate and period 
orbit_rate=sqrt(earth_gravity_constant/radius‘’3); 
period=2*pi/orbit_rate; 


% DCM relating perifocal PQW frame to Earth-centered ECI frame 
eci_from_pqw=[cos(ra) -sin(ra)*cos(inc) sin(ra)*sin(inc) 
sin(ra) cos(ra)*cos(inc) -cos(ra)*sin(inc) 
0 sin(inc) cos(inc)]; 


% initial Greenwich sidereal time 
greenwich_time=0*d2r; 


% satellite inertia matrix and its inverse 
satellite inertia=[3000 %***** UPDATE DATA ITEM 
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0310 
0 0 2]; 
inv_satellite_inertia=inv(satellite inertia); 


% satellite torque coil capability 
m_max=30; MW***EE UPDATE DATA ITEM 


dipole_off=input('Disable mag torquers? (1=Yes,0=No)'); 


% magnetic b-dot control feedback gains 
k_mag=100000000; 


% initial satellite state 

true_anomaly=0*d2r; 

“%rates=[0.01;0.01;0.01]; 

rates= [-0.1; 0.1; 0.01]; % NPS PDR sim results used these rates 

Y%rates= [ 0.0; 0.001083; 0.0]; % orbital rate @ 600km altitude on Y axis 

“%quats=[0;0;0;1]; 

quats= [0.62721137512625; 0.32650557562198; % ran routine initquat.m 
0.62721137512625; 0.32650557562198]; % to get this 


state=[rates;quats]; 
bfield_xyz_old=[0;0;0]; 


% begin time loop 

dt= 2; M***EE UPDATE DATA ITEM 
%t_final=3600*.5; 

t_final= input(Input sim duration (seconds): '); % total sim time 


for t=0:dt:t_ final 
% DCM relating spacecraft XYZ frame to ECI frame 


xyz_from_eci=quats_to_dcm(quats); 


% DCM relating rotating orbit frame 


% (x along velocity, z along zenith) to PQW frame 
orb_from_pqw=[-sin(true_anomaly) cos(true_anomaly) 0 
1 


cos(true_anomaly) sin(true_ anomaly) 0]; 


% DCM relating rotating orbit frame to ECI frame 
orb_from_eci=orb_from_pqw*eci_from_pqw’; 


% DCM relating spacecraft XYZ frame to rotating orbit frame 
xyz_from_orb=xyz_from_eci*orb_from_eci'; 


% 3(yaw)-2(pitch)-1(roll) Euler angle sequence describing 
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% spacecraft frame relative to orbit frame 
yaw=atan(xyz_from_orb(1,2)/xyz_from_orb(1,1)); 
pitch=asin(-xyz_from_orb(1,3)); 
roll=atan(xyz_from_orb(2,3)/xyz_from_orb(3,3)); 


% position of satellite in PQW frame 
r_pqw=radius*[cos(true_anomaly);sin(true_anomaly);0]; 


% position of satellite in ECI frame 
r_eci=eci_from_pqw*r_pqw; 





% position of satellite in ECEF frame 
ecef from_eci=[cos(greenwich_time) sin(greenwich_time) 0 
-sin(greenwich_ time) cos(greenwich_time) 0 
0 0 1]; 
r_ecef=ecef_from_eci*r_eci; 


% gravity gradient torque 

r_xyz=xyz_from_eci*r_eci; 
r_xyz_unit=r_xyz/sqrt(dot(r_xyz,r_xyz)); 
u_ge=3*orbit_rate*2*skew(r_xyz_unit)*satellite_inertia*r xyz_unit; 


% angle of spacecraft z-axis relative to local zenith direction (libration angle) 
theta=acos(r_xyz_unit(3)); 


% local magnetic field in ECEF frame 
bfield_ecef=bfield(r_ecef,2); 


% local magnetic field in XYZ frame 
bfield_xyz=xyz_ from eci*ecef_from_eci'*bfield_ecef; 


% time derivative of local magnetic field in XYZ frame 
bdot=(bfield_xyz-bfield_xyz_old)/dt; 
bfield_xyz_old=bfield_xyz; 


% commanded spacecraft dipole 


if dipole_off==1, 
dipole(1:3)=0; 
else 


if t < 3600*30 % may want to remove this time-out 
dipole(1)=-k_mag*bdot(1); 
dipole(2)=-k_mag*bdot(2); 
dipole(3)=-k_mag*bdot(3); 

else 
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dipole(1:3)=0; 
end 
end 


% dipole saturation 
for 1=1:3 
if abs(dipole(i)) > m_max 
dipole(i)=m_max*sign(dipole(i)); 
end 
end 


% corresponding magnetic torque 
u_mag=cross(dipole,bfield_xyz)'; 


% store critical outputs before updating 

count=count+1; 

if count == store_count 
k=k+1; 
count=0; 
time(k)=t; 
w1(k)=rates(1); 
w2(k)=rates(2); 
w3(k)=rates(3); 
m1(k)=dipole(1); 
m2(k)=dipole(2); 
m3(k)=dipole(3); 
b1l(k)=bfield_xyz(1); 
b2(k)=bfield_xyz(2); 
b3(k)=bfield_xyz(3); 
ugg1(k)=u_gg(1); 
ugg2(k)=u_gg(2); 
ugg3(k)=u_gg(3); 
umag1(k)=u_mag(1); 
umag2(k)=u_mag(2); 
umag3(k)=u_mag(3); 
yaw_angle(k)=yaw; 
pitch_angle(k)=pitch; 
roll_angle(k)=roll; 
libration(k)=theta; 

end 


% ****E PROPAGATE ATTITUDE AND ORBIT STATES ***** 


% updated rotational state of satellite relative to ECI frame 
u_ext=u_ggtu_mag; 
state=satellite_dynamics(state,t,u_ext); 
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rates=state(1:3); 
quats=state(4:7); 
quats=qnorm(quats); 


% updated greenwich time 
greenwich time=greenwich_time+earth_rate*dt; 


% updated true anomaly 
true_anomaly=true_anomalytorbit_rate*dt; 
end 


plotresult 


function state=satellite_dynamics(state_old,t,u_ext) 
global satellite inertia inv_satellite inertia dt 


rates_temp=state_old(1:3); 

q_temp=state_old(4:7); 

rates_dot=inv_satellite_inertia* 
(u_ext-skew(rates_temp)*satellite_inertia*rates_temp); 

qmatrix=[q_temp(4) -q_temp(3) q_temp(2) q_temp(1); 
q_temp(3) q_temp(4) -q_temp(1) q_temp(2); 
-q_temp(2) q_temp(1) q_temp(4) q_temp(3); 
-q_temp(1) -q_temp(2) -q_temp(3) q_temp(4)]; 


quats_dot=0.5*qmatrix*[rates_temp;0]; 
fl=dt*[rates_dot;quats_ dot]; 


rates temp=state_old(1:3)+f1(1:3)/2; 
q_temp=state_old(4:7)+f1(4:7)/2; 
rates_dot=inv_satellite_inertia* 
(u_ext-skew(rates_temp)*satellite_inertia*rates_temp); 
qmatrix=[q_temp(4) -q_temp(3) q_temp(2) q_temp(1); 
q_temp(3) q_temp(4) -q_temp(1) q_temp(2), 
-q_temp(2) q_temp(1) q_temp(4) q_temp(3); 
-q_temp(1) -q_temp(2) -q_temp(3) q_temp(4)]; 


quats_dot=0.5*qmatrix*[rates_temp;0]; 
f2=dt*[rates_dot;quats_ dot]; 


rates _temp=state_old(1:3)+f2(1:3)/2; 

q_temp=state_old(4:7)+f2(4:7)/2; 

rates_dot=inv_satellite_inertia* 
(u_ext-skew/(rates_temp)*satellite_inertia*rates_temp); 

qmatrix=[q_temp(4) -q_temp(3) q_temp(2) q_temp(1); 
q_temp(3) q_temp(4) -q_temp(1) q_temp(2); 
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-q_temp(2) q_temp(1) q_temp(4) q_temp(3); 

-q_temp(1) -q_temp(2) -q_temp(3) q_temp(4)]; 
quats_dot=0.5*qmatrix*[rates_temp;0]; 
f3=dt*[rates_dot;quats_dot]; 


rates _temp=state_old(1:3)+f3(1:3); 

q_temp=state_old(4:7)+f3(4:7); 

rates _dot=inv_satellite_inertia*(u_ext-skew(rates_temp)* 

satellite inertia*rates_ temp); 

qmatrix=[q_temp(4) -q_temp(3) q_temp(2) q_temp(1); 
q_temp(3) q_temp(4) -q_temp(1) q_temp(2); 
-q_temp(2) q_temp(1) q_temp(4) q_temp(3); 
-q_temp(1) -q_temp(2) -q_temp(3) q_temp(4)]; 

quats_dot=0.5*qmatrix*[rates_temp;0]; 

f4=dt*[rates_dot;quats_dot]; 


state=state_old+(f1+2*f2+2*f3+f4)/6; 





function bfield_ecef=bfield(r_ecef,order) 
global earth_radius 


if order > 8 

order=8; 

‘Desired order exceeds 8 --- model order reset to 8' 
end 


% current orbit position in Earth-fixed REF frame (up, south, east) 
r=sqrt(r_ecef(1)42+1_ecef(2)2+1_ecef(3)2); 
delta=asin(r_ecef(3)/r); 
cdelta=cos(delta); 
sdelta=sin(delta); 
theta=pi/2-delta; 
stheta=sin(theta); 
ctheta=cos(theta); 
phi=acos(r_ecef(1)/(r*cdelta)); 
sphi=r_ecef(2)/(r*cdelta); 
cphi=r_ecef(1)/(r*cdelta); 
if sphi<0 

phi=-phi; 
end 


ref_from_ecef=[cdelta*cphi cdelta*sphi sdelta 
sdelta*cphi sdelta*sphi -cdelta 
-sphi cphi 0]; 
% S factors 
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s0(1)=1; 
s(1,1)=1; 
if order > 1 
for n=2:order 
s0(n)=s0(n-1)*(2*n-1)/n; 
s(n, 1)=s0(n)*sqrt(2*n/(n+1)); 
for m=2:n 
s(n,m)=s(n,m-1)*sqrt((n-m+1)/(n+m)); 
end 
end 
end 


% Legendre functions and their derivatives wrt theta 
p0(1)=ctheta; 
p(1,1)=stheta; 
dp0(1)=-stheta; 
dp(1,1)=ctheta; 
if order > 1 
for n=2:order 
if n== 
p0(2)=ctheta’2-1/3; 
dp0(2)=-2*stheta*ctheta; 
else 
p0(n)=p0(n-1)*ctheta-(n-1)*2/((2*n-1)*(2*n-3))*p0(n-2); 
dp0(n)=dp0(n-1)*ctheta-p0(n-1)* 
stheta-(n-1)*2/((2*n-1)*(2*n-3))*dp0(n-2); 
end 
p(n,n)=p(n-1,n-1)*stheta; 
dp(n,n)=dp(n-1,n-1)*stheta+p(n-1,n-1)*ctheta; 
for m=1:(n-1) 
if m==(n-1) 
p(n,m)=p(n-1,m)*ctheta; 
dp(n,m)=dp(n-1,m)*ctheta-p(n-1,m)*stheta; 
else 
p(n,m)=p(n-1,m)*ctheta-((n-1)*2-m%2)/((2*n-1)*(2*n-3))*p(n-2,m); 
dp(n,m)=dp(n-1,m)*ctheta-p(n-1,m)* 
stheta-((n-1)2-m%2)/((2*n-1)*(2*n-3))*dp(n-2,m); 
end 
end 
end 
end 


% look-up table for Gaussian coefficients 
coeffs=gauss_coeffs; 

g=coeffs(1:44); 

h=coeffs(45:88); 
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% first order, zeroth degree field components 
br=2*(earth_radius/r)*3*s0(1)*g(1)*p0(1); 
bt=-(earth_radius/r)*3*s0(1)*g(1)*dp0(1); 
bp=0; 


% higher orders, zeroth degree field components 
if order > 1 
Jale 
for n=2:order 
jt; 
br=br+(n+1)*(earth_radius/r)(n+2)*sO0(n)*g(q)*p0(n); 
bt=bt-(earth_radius/r)*(n+2)*sO(n)*g(j)*dp0(n); 
end 
end 


% all orders, higher degree field components 
j=l; 
for n=1:order 
jal 
for m=1:n 
br=br+(n+1)*(earth_radius/r)(n+2)*s(n,m)*(g(j)* 
cos(m*phi)+h(j)*sin(m*phi))*p(n,m); 
bt=bt-(earth_radius/r)(n+2)*s(n,m)*(g(j)*cos(m*phi)+h(Gj)* 
sin(m*phi))*dp(n,m); 
bp=bp-(1/stheta)*(earth_radius/r)(n+2)*m* 
s(n,m)*(-g(j)*sin(m*phi)+h(j)*cos(m*phi))*p(,m); 
i ela 
end 
end 


% B-field in REF frame 
bfield_ref=[br;bt;bp]; 

% B-field in ECEF frame 
bfield_ecef=ref_from_ecef'*bfield_ref; 


function coeffs=gauss_coeffs 


g=[-29682 -1789 -2197 3074 1685 1329 -2268 1249 769 941 782 291 —421 
116 -210 352 237 -122 -167 -26 66 64 65 -172 2 17 -94 78 -67 
12948 10 -2 244 -1 -9 -1445 0 -7]*1.0e-09; 
h=[0 5318 0 -2356 -425 0 -263 302 -406 0 262 -232 98 -301 0 44 157-152 
-64 99 0 -16 77 67 -57 4 28 0 -77 -25 3 22 16 -23 -3 0 12 -207 —21 
12 10 -17 -10]*1.0e-09; 


coeffs=[g h]; 
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function qn=qnorm(q) 


mag=sqrt(dot(q,q)); 
qn=[q(1);q(2);q(3);q(4) ]/mag; 


function dem=quats_ to _dcm(q) 
dem=[q(4)*2+q(1)%2-q(2)%2-q(3)%2 2*(q(1)*q(2)+q(3)*q(4)) 
2*(q(1)*q(3)-q(2)*q(4)); 


2*(q(1)*q(2)-q(3)*q(4)) q(4)*2-q(1)°2+q(2)°2-q(3)"2 
2*(q(2)*q(3)+q(1)*q(4)); 
2*(q(1)*q(3)+q(2)*q(4)) 2*(q(2)*q(3)-q(1)*q(4)) 


q(4)*2-q(1)*2-q(2)2+q(3)*2]; 


function skew_matrix=skew(vector) 


skew_matrix=[0 -vector(3) vector(2); 
vector(3) 0 -vector(1); 
-vector(2) vector(1) Os 


“%plotresult.m 
r2d= 180/pi; 


figure 

subplot(311),plot(time,roll_angle*r2d),title('Attitude Error') 
ylabel('Roll(deg)'),grid 

subplot(312),plot(time,pitch_angle*r2d), ylabel('Pitch(deg)'), grid 
subplot(313),plot(time,yaw_angle*r2d), ylabel("Yaw(deg)'), grid 
xlabel(‘time(sec)') 


figure 

subplot(311),plot(time,w 1 *r2d),title('Body Rates’), ylabel('wx(deg/sec)'),grid 
subplot(3 12),plot(time,w2*r2d), ylabel('wy(deg/sec)'), ylabel('wy(deg/sec)'), grid 
subplot(3 13), plot(time,w3*r2d), ylabel('wz(deg/sec)'), ylabel('wz(deg/sec)'), grid 
xlabel(‘time(sec)') 


NPS (LABISMOIDATA.M) 
% Lab 15MOlIdata.m written by Prof Barry Leonard 





Ix=5;ly=5.1;Iz=2;kax=.028;kay=.015;kaz=.032*4; % save for b8 
% steady state mag act 





Ixy=0;Ixz=0;lyz=0; 

beta=0*pi/180; 

Imoi=[Ix -Ixy -Ixz;-Ixy ly -lyz;-Ixz -lyz Iz]; 
linv=inv(Imoi); 


ral 








pho=-.1;tho=.1;pso=-.1; 

phdo=-0. 1;thdo=0.1;psdo=0.01; 
s1=sin(pho/2);s2=sin(tho/2);s3=sin(pso/2); 
cl=cos(pho/2);c2=cos(tho/2);c3=cos(pso/2); 
elo=sl*c2*c3-cl*s2*s3;e20=cl*s2*c3+s1*c2*s3; %Wie pg.321 
e30=c1 *c2*s3-s1*s2*c3;e40=c 1 *c2*c3+s1*s2*s3; 
S1=sin(pho);S2=sin(tho);S3=sin(pso); 
C1=cos(pho);C2=cos(tho);C3=cos(pso); 





wxo=phdo-psdo*S2-wo*S3*C2; 
wyo=thdo*C1+psdo*C2*S 1-wo*(C3*C1+S83*S2*S 1); 
wzo=psdo*C2*C1-thdo*S 1-wo*(S3*S2*C1-C3*S1); 
eo=[elo e20 e30 e40];qo=e0; 
Ho=Imoi*[wxo wyo wzo]'; 


sat=30; 
g1=0.75;g2=0.82;¢3=0.43; % varies with inclination 

% (from torque ave. with high sat.) 
ao=Ix-ly+1z;a1=ly-Iz;a2=Ix-Iz;a3=ly-Ix; 


A=[000100;000010;000001; 
-4*wo%2*al/Ix 0 0 0 0 wo*ao/Ix; 
0 -3* wo’2*a2/ly 00 0 0; 
0 0 -wo%2*a3/Iz -wo*ao/Iz 0 0]; 
B=Kme*[0 0 0;0 0 0;0 0 0;kax/Ix 0 0;0 kay/Iy 0;0 0 kaz/Iz]; 


Qx=110*eye(6); Ru=1*[.1 0 0;0 .095 0;0 0 .008];%for 5 5.1 2 with b8 
[K,S,e]-lqr(A,B,Qx,Ru); 


Aaa=A(1:3,1:3);Aab=A(1:3,4:6);Aba=A (4:6, 1:3); Abb=A(4:6,4:6); 
Ba=B(1:3,1:3);Bb=B(4:6,1:3);Ka=K(1:3,1:3);Kb=K(1:3,4:6); 
Ga=G(1:3,1:3);Gb=G(4:6, 1:3); 

Lr=1*diag([1 .5 .6]); %for 5 5.1 2 with b8& 

Lr=0.09* diag([1 .5 .6]); 


Ar=Abb-Lr* Aab-(Bb-Lr*Ba)*Kb; 
Br=Ar*Lr+Aba-Lr* Aaa-(Bb-Lr*Ba)*Ka; 
Cr=-Kb;Dr=-Ka-Kb* Lr; 

Kp=Ka; Kd=Kb; 


s1=sin(pho/2);s2=sin(tho/2);s3=sin(pso/2); 
cl=cos(pho/2);c2=cos(tho/2);c3=cos(pso/2); 
elo=sl*c2*c3-cl*s2*s3;e20=c 1 *s2*c3+s1*c2*s3;%Wie pg.321 
e30=c1 *c2*s3-s1*s2*c3;e40=c 1 *c2*c3+81*s2*s3; 
S1=sin(pho);S2=sin(tho);S3=sin(pso); 





a2 


C1=cos(pho);C2=cos(tho);C3=cos(pso); 


wxo=phdo-psdo*S2-wo*S3*C2; 
wyo=thdo*C1+psdo*C2*S 1-wo*(C3*C1+S83*S2*S 1); 
wzo=psdo*C2*C1-thdo*S 1-wo*(S3*S2*C1-C3*S1); 
eo=[elo e20 e30 e40];qo=e0; 
Ho=Imoi*[wxo wyo wzo]'; 


NPS (LAB16DATA.M) 
“%Lab16Data written by Prof Barry Leonard 


clear 


Re=6371.2e3;mu=398601.2e9;we=7.292 1e-5;epsilon=11.398*pi/180; 
Altitude=[450 500 550 600 650 700 750 800]*1e3; “look up table data for aero 
Density=[36.1 18 9.25 4.89 2.64 1.47 .837 .439]*le-13; % look up table 

% data for aero 


h=600e3; % spacecraft altitude (km) 
incln=(35)*pi/180; % inclination (radians) 
beta=0*pi/180; % angle between outward normal vector and sun vector 


nuo=-115*pi/180; 

alphao=0*pi/180; 

Lm=-70.454;Lgo=0;u0=(Lgo+Lm+90)*pi/180; 
wn=0;wo=sqrt(mu/(Reth)‘3); 
V=wo*(Reth);rho=asin(Re/(h+Re));Kme=7.9638e15/(Reth)43; 
Cd=2.5;Kaero=0.5*Cd*V%2; % Cd = drag coefficient 
psun=4.5e-6;Psolar=2*psun; 

Area=[0.2674 0.2674 .1927]; dL=[.002 .002 .008]; 


Ixy=0;Ixz=0;lyz=0; 
Ix=5;Iy=5.1;Iz=2;kax=1;kay=0.1*kax;kaz=kax* 1; 








Imoi=[Ix -Ixy -Ixz; % The Inertia tensor or Inertia Matrix 

-Ixy ly -lyz; 

-Ixz -lyz Iz]; 
linv=inv(Imoi); % Inertia Matrix Inverse 
sat=30; 


g1=0.75;g2=0.82;¢3=0.43; % varies with inclination 
% (from torque ave. with high sat.) 


ao=Ix-ly+1z;a1=ly-Iz;a2=Ix-Iz;a3=ly-Ix; 

A=[000100;0000 1 0;0 0000 1;-4*wo’2*al/Ix 0 0 0 0 wo*ao/Ix;... 
0 -3*wo%2*a2/Ty 0 0 0 0;0 0 -wo’2*a3/Iz -wo*ao/Iz 0 0]; 

B=Kme*[0 0 0;0 0 0;0 0 0;kax/Ix 0 0;0 kay/Ty 0;0 0 kaz/Iz]; 

G=[0 0 0;0 0 0;0 0 0;1/Ix 0 0;0 1/Iy 0;0 0 1/Iz]; 

Qx=10*eye(6); Ru=[.1 0 0;0 .1 0;0 0 .01]; 
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[K,S,e]=lqr(A,B,Qx,Ru); 


Aaa=A(1:3,1:3);Aab=A(1:3,4:6);Aba=A (4:6, 1:3); Abb=A(4:6,4:6); 
Ba=B(1:3,1:3);Bb=B(4:6,1:3);Ka=K(1:3,1:3);Kb=K(1:3,4:6); 
Ga=G(1:3,1:3);Gb=G(4:6, 1:3); 

Lr=0.2*eye(3);Ar=Abb-Lr* Aab-(Bb-Lr*Ba)*Kb; 
Br=Ar*Lr+Aba-Lr* Aaa-(Bb-Lr*Ba)*Ka; 

Cr=-Kb;Dr=-Ka-Kb* Lr; 








Kp=Ka; Kd=Kb; 
pho=-.1;tho=.1;pso=-.01; % Euler angles 
phdo=.1;thdo=-0.01;psdo=0.01; % Euler rates 


s1=sin(pho/2);s2=sin(tho/2);s3=sin(pso/2); 
cl=cos(pho/2);c2=cos(tho/2);c3=cos(pso/2); 


elo=sl*c2*c3-cl*s2*s3;e20=cl*s2*c3+s1*c2*s3; % Wie pg.321 
e30=c1 *c2*s3-s1*s2*c3;e40=c 1 *c2*c3+s1*s2*s3; 





S1=sin(pho);S2=sin(tho);S3=sin(pso); 
C1=cos(pho);C2=cos(tho);C3=cos(pso); 


wxo=phdo-psdo*S2-wo*S3*C2; % WX 
wyo=thdo*C1+psdo*C2*S 1-wo*(C3*C1+S3*S2*S1);  % wy 
wzo=psdo*C2*C1-thdo*S 1-wo*(S3*S2*C1-C3*S 1); % WZ 


eo=[elo e20 e30 e40]; 


qo=eo; 

Ho=Imoi*[wxo wyo wzo]'; % [Ho] = Angular Momentum [3x1] 

“%SphericalHarmMag3 

r=Reth;a=Re; % r =radius from sat to center of earth 
% Re = a= Earth radius = 6371.2e3 km 

kme=Kme* 1e9; %conversion from nT.m‘%3 to T.m43 


dth=5;tmin=0;tmax=180;Jmax=(tmax-tmin)/dth; 
dphi=10;phimin=-180;phimax=180;Imax=(phimax-phimin)/dphi; 


gge=[-30186 -2036 00 0000 0;-1898 2997 15510000003... 
1299 -2144 1296 805 00 00 0;951 807 462 -393 235 00 0 03... 
-204 368 275 -20 -161 -38 0 0 0;46 57 15 -210 -1 -8 -114 0 03... 
66 -57 -7 7 -22 -9 11 -8 0;11 13 3 -12 -46 -29 1]; 


hhh=[0 5735 0 00 0 0 0 0;0 -2124 -3700000 0... 
0 -361 249 -253 0000 0;0 148 -264 37 -307 00 0 03... 
0 39 142 -147 -99 74 0 0 0;0 -23 102 88 -43 -9 -4 0 03... 
0 -68 -24 -4 11 27 -17 -14 0;0 4 -15 2 -19 1 18 -6 -19]; 
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gg=[-29615 -1728 00 0 0 0 0 0;-2267 3072 167200000 03... 
1341 -2290 1253 715 000 0 0;935 787 251 -405 110000 03... 
-217 351 222 -131 -169 -12 0 0 0;72 68 74 -161 -5 17 -91 0 03... 
79 -74 033 97 8 -2 0325 6 -9 -8 -17 97 -8 -7]; 


hh=[0 5186 000000 0;0 -2478 -45800000 03... 
0 -227 296 -492 00 0 0 0;0 272 -232 119 -304 00 0 03... 
0 44 172 -134 -40 107 0 0 0;0 -17 64 65 -61 1 44.0 03... 
0 -65 -24 6 24 15 -25 -6 0;0 12 -22 8 -21 15 9 -16 -3]; 


for J=1:Jmax+1; 
t=(J-1)*dth*pi/180; th(J)=(J-1)*dth; 
Ct=cos(t);St=sin(t); 


for I=1:Imax+1; 
p=(phimin+(I-1)*dphi)*pi/180; phi(1)=phimin+(1-1)*dphi; 
Cp=cos(p);Sp=sin(p); 


for i=1:8 
for j=1:i+1 
Km(i,j)=((i-1)42-(G-1)42)/(2*1-1)/(2*1-3); 
if i==1; Km=0;end 
end 
end 


Km=tril(Km, 1); 
for i=1:8 
for j=1:i+1 
P(1,1)=Ct; 


if i==1 & j==2;P(i,j)=St* 1send 
if i> 1 & j==i+1;P(,j)=St*P(i-1,j-1);end 
if i==2;Pk=1;end; 
if i>2;Pk=P(i-2,j);end 
if i>1 & j~=it+1;PC,j)=Ct*P-1,j)-Km(i,j)*Pk;end 
end 
end 


P=tril(P,1); 
for i=1:8 
for j=1:i+1 
dp(1,1)=-St; 
if i==1 & j==2;dp(i.j)=St*0+Ct* 1 end 
if i>1 & j==i+1;dp(ij)=St* dp(i-1,j-1)+Ct*P(i-1,j-1);end 
if i==2;dpk=0;end; 
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if i>2;dpk=dp(i-2,j);end 
if i> 1 & j~=i+1;dp(,j)=Ct* dp(i-1,j)-St*P(-1,))-Km(i,j)*dpk;end 
end 

end 


dp=tril(dp,1); 
for i=1:8 
for j=1:9 
dm=0;if j==2;dm=1;end 
if j==1 & i==1; S=l;end 
if j==1 & 1-=1;S(,jJ)=SG-1,j)*(2*i- 1D /i;end 
if j>1;SG,j)=SGJ-1)*sqrt(-G-1)+1)*(dm+1)/(i+G-1)));end 
end 
end 














S=tril(S,1); 
G=S."*gg; 
H=S.*hh; 
for i=1:8 
j=l:i+1;m=j-1; 
F(i,j)=G(i,j).*cos(m*p)+H(i,j).*sin(m*p); 
end 


for 1=1:8 
j=l:i+1;m=j-1; 
f(ij)=m.*(-G(i,j).*sin(m*p)+H(i,j).*cos(m*p)); 
end 


if St==0; St=0.001; ;end 

for i=1:8 
aal (i)=(i+1)*(a/r)(it2); 
aa2(1)=-aa 1 (i)/(i+1); 
aa3(i)=aa2(i)/St; 

end 


FPS=sum((F.*P)')' ; Bl=aal'.*FPS; Br8=sum(B1); Brl=B1(1); 
FdpS=sum((F.*dp)')'; B2=aa2'.*FdpS; Bt8=sum(B2); Bt2=B2(1); 
fPS=sum((f.*P)')' ; B3=aa3'.*fPS; Bp8=sum(B3); Bp3=B3(1); 


br8(1,J)=Br8/kme; br1(,J)=Br1/kme; 
bt8(1,J)=Bt8/kme;  bt1(1,J)=Bt2/kme; 
bp8(1,J=Bp8/kme; bp1(1,J)=Bp3/kme; 
b1(L,J=(Brl.42+Bt2..2+Bp3.2).%0.5/kme; 
b8(L,J=(Br8.42+Bt8..2+Bp8.%2).%0.5/kme; 
end 

end 
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Inclination=[0.01 10 20 30 40 50 60 70 80 90 100 110 120]*pi/180; 
gx=[.967 .955 .922 .876 .826 .78 .739 .709 .691 .686 .691 .709 .739]; 
gy=[0.0001 .0678 .256 .46 .632 .762 .857 .923 .965 .981 .965 .923 .857]; 
gz=[.804 .781 .711 .614 522 .446 .39 .353 .335 .333 .335 .353 .39]; 


NPS (LAB17.MDL) 


See Figure Al. 
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NPS (NPSATIACSDATA.M) 
%NPSAT1ACSData written by Prof Barry Leonard 


clear 

tic 

Re=6371.2e3;mu=398601.2e9;we=7.2921e-5;%earth radius, gravity and spin rate 
Altitude=[450 500 550 600 650 700 750 800]*1e3;%look up table data for aero 
Density=[36.1 18 9.25 4.89 2.64 1.47 .837 .439]*1le-13;%density table data 


h=600e3;incln=(35)*pi/180;beta=15*pi/180;%altitude,inclination, solar angle 
nuo=-115*pi/180;alphao=0*pi/180;%initial subsolar point and true anomaly 
Lgo=0;%initial position of Greenwich Meridian WRT RAAN 

kpre=0; “onodal precession constant assumed to be zero here 
wn=kpre*(Re/(Ret+h))*3.5*cos(incln);“%onodal precession (zero eccentricity) 
wo=sqrt(mu/(Re+h)*3);V=wo*(Reth); “orbital angular and linear velocity 
rho=asin(Re/(h+Re)); Y%earth angular radius 


Cd=2.5;psun=4.5e-6;%odrag coefficient and solar pressure constant~N/m‘2 
Kaero=0.5*Cd*V“2;Psolar=2*psun; % constants for aero and solar torque calc. 


Area=[0.2674 0.2674 .1927]; Yoprojected area~m”’2 in body x,y,z directions 
dL=[.002 .002 .008]; Yopredicted (cp-cm)~m in body x,y,z directions 








Ix=5;ly=5.1;Iz=2;Ixy=0;Ixz=0;Iyz=0; % moments of inertia (MOI)~kg.m“2 
Imoi=[Ix -Ixy -[xz;-Ixy ly -lyz;-Ixz -lyz Iz]; %MOI matrix 

Tinv=inv(Imoi); %MOI matrix inverse 
ao=Ix-ly+Iz;al=ly-Iz;a2=Ix-Iz;a3=Iy-Ix;%MOI combination constants 








pho=-.1;tho=.1;pso=-0.1; Yoinitial Euler angles (phi,theta,psi) ~r 
phdo=-0.1;thdo=0.1;psdo=0.01; % initial Euler angle rates ~r/s 


% calculation of initial quaternion (qo) and angular momentum (Ho): 
s1=sin(pho/2);s2=sin(tho/2);s3=sin(pso/2);c 1=cos(pho/2);c2=cos(tho/2); 
c3=cos(pso/2); qlo=s1*c2*c3-cl*s2*s3;q20=cl*s2*c3+s1*c2*s3;% Wie pg.321 
q30=c1*c2*s3-s1*s2*c3;q40=cl*c2*c3+s1*s2*s3; 
S1=sin(pho);S2=sin(tho);S3=sin(pso);C 1=cos(pho);C2=cos(tho);C3=cos(pso); 
wxo=phdo-psdo*S2-wo*§S3*C2; 
wyo=thdo*C1+psdo*C2*S 1-wo*(C3*C1+S83*S2*S 1); 
wzo=psdo*C2*C1-thdo*S 1-wo*(S3*S2*C1-C3*S 1); 

qo=[qlo q20 q3o0 q4o]; 

Ho=Imoi*[wxo wyo wzo]'; 


“Calculation of spherical harmonic magnetic field model (Wertz pp779-783) 
t=Reth;a=Re; “definitions to match Wertz pg 780 
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% IGRF Epoch 2000 Guassian coefficients ~n.T : 
gg=[-29615 -1728 0000 0 0 05-2267 3072 16720000003... 
1341 -2290 1253 715 0000 0;935 787 251 -405 110000 03... 
-217 351 222 -131 -169 -12 0 0 0;72 68 74 -161 -5 17 -91 0 03... 
79 -74 0 33 97 8 -2 0;25 6 -9 -8 -17 9 7 -8 -7]; 
hh=[0 5186 0 0000 0 0;0 -2478 -45800000 0... 
0 -227 296 -492 00 0 0 0;0 272 -232 119 -30400 0 03... 
0 44 172 -134 -40 107 0 0 0;0 -17 64 65 -61 1 44.003... 
0 -65 -24 6 24 15 -25 -6 0;0 12 -22 8 -21 15 9 -16 -3]; 
Kme=(a/r)3*sqrt(gg(1,1)42+gg(1,2)2+hh(1,2)2)* le-9;Y%odipole strength~Tesla 


ka=[.02863 .01534 .1309];sat=30;%actuator gains and saturation 
%ka=[.03 .015 .13]; 
%Calc of state space A,B,G matrices follows (xdot=Ax+Bu+Gw,w=disturbance) 
A=[000100;0000 1 0;0 0000 1;-4*wo’2*al/Ix 0 0 0 0 wo*ao/Ix;... 
0 -3*wo%2*a2/ly 0 0 0 0;0 0 -wo’%2*a3/Iz -wo*ao/Iz 0 0]; 
B=Kme*[0 0 0;0 0 0;0 0 0;ka(1)/Ix 0 0;0 ka(2)/Iy 0;0 0 ka(3)/Iz]; 
G=[0 0 0;0 0 0;0 0 0;1/Ix 0 0;0 1/Iy 0;0 0 1/Iz]; 


Qx=110*eye(6); Ru=1*[.1 0 0;0 .095 0;0 0 .008]; %LQR weighting matrices 
[K,S,e]=lqr(A,B,Qx,Ru); %LQR gain calculation 


“partitioning of A,B,G matrices required for reduced order estimator: 
Aaa=A(1:3,1:3);Aab=A(1:3,4:6);Aba=A (4:6, 1:3); Abb=A(4:6,4:6); 
Ba=B(1:3,1:3);Bb=B(4:6, 1:3);Ka=K(1:3,1:3);Kb=K(1:3,4:6); 
Ga=G(1:3,1:3);Gb=G(4:6, 1:3); 

Lr=0.09* diag([1 .5 .6]); %Estimator gain found by simulation 


del_t=5;tmin=0;tmax=180;%colatitude increment, min, max (theta) 
Jmax=(tmax-tmin)/del_t; % max value of J index 
del_p=10;pmin=-180;pmax=180; “longitude increment, min, max (phi) 
Imax=(pmax-pmin)/del_p; %omax value of I index 


for J=1:Jmax+1; 
t=(J-1)*del_ t*pi/180;Ct=cos(t);St=sin(t); Yosine & cosine of theta 
th(J)=(J-1)*del_t; Yotheta vector (look up tables Mag Fld Model) 
for I=1:Imax+1; 
p=(pmin-+(I-1)*del_p)*pi/180;Cp=cos(p);Sp=sin(p);%osin & cos phi 
phi(I)=pmin+(I-1)*del_p;%phi vector (look up tables Mag Fld Model) 
% Calculation of Legendre function constant (Km) follows 
for i=1:8 
for j=1:i+1 
Km(i,j)=((i-1)42-(Gj-1)42)/(2*1-1)/(2*1-3); 
if i==1; Km=0;end 
end 
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end 
Km=tril(Km, 1); 


%Calculation of Legendre function (P) follows 
for i=1:8 
for j=1:i+1 
P(1,1)=Ct; 
if i==1 & j==2;P(i,j)=St* 1;end 
if i>1 & j==i+1;P(,j)=St*P(i-1,j-1);end 


if i==2;Pk=1;end; 
if 1>2;Pk=P(i-2,j);end 
if i>] & j~=i+1;PC,j)=Ct*P-1,j))-Km(i,j)*Pk;end 
end 
end 
P=tril(P,1); 


% Calculation of partial of Legendre function WRT theta (dp) follows 
for 1=1:8 
for j=1:i+1 
dp(1,1)=-St; 
if i==1 & j==2;dp(i,j)=St*0+Ct* 1 end 
if i> 1 & j==i+1;dp(i,j)=St* dp(-1,j-1)+Ct*P(-1,j-1);end 
if i==2;dpk=0;end; 
if i>2;dpk=dp(i-2,j);end 
if i> 1 & j~=i1+1;dp(,j)=Ct*dp(i-1,j)-St*P(-1,))-Km(i,j)*dpk;end 
end 
end 
dp=tril(dp,1); 





% Calculation of Gausian coefficient norm. factor (S) follows 
for i=1:8 
for j=1:9 
dm=0;if j==2;dm=1;end 
if j==1 & 1==1; S=1;end 
if j==1 & 1-=1;S(,jJ)=SG-1,))*(2*i-1/i;end 
if > 1;SG,)=S(GJ-1)*sqrt(a-G-1)+1)*(dm+1)/(i+Q-1)));end 
end 
end 
S=tril(S,1); 














%Calculation of normalized Gaussian coefficients ~ n.T : 
Geg=S.* gg; Hh=S.*hh; 


for 1=1:8 
jul:i+1;m=j-1; 
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F(i,j)=Gg(i,j).*cos(m*p)+Hhi(,j).*sin(m*p); Yodefinition 
end 
for i=1:8 
jel:it+1j;m=j-1; 
f(i,j)=m.*(-Gg(i,j).*sin(m*p)+Hh(i,j).*cos(m*p));% definition 
end 


if St==0; St=0.001; ;end “singularity avoidance in following loop 
% Calculation of altitude dependent coefficients follows 
for i=1:8 

aal (i)=(i+1)*(a/r)(it+2); aa2(i)=-aal (1)/(i+ 1); aa3(i)=aa2(i)/St; 
end 


% Calculation of radial,theta and phi magnetic field components ~Tesla 
FPS=sum((F.*P)')' ; Bl=aal'.*FPS; Br8=sum(B1)*1e-9; 
FdpS=sum((F.*dp)')'; B2=aa2'.*FdpS; Bt8=sum(B2)*1e-9; 
fPS=sum((f.*P)')' ; B3=aa3'.*fPS; Bp8=sum(B3)*le-9; 


“Calculation of normalized field components (input to lookup tables) 
br8(1,J)=Br8/Kme; 

bt8(,J)=Bt8/Kme; 

bp8(1,J)=Bp8/Kme; 


% Calculation of field vector normalized magnitude: 
b8(1,J=(Br8.42+Bt8.*2+Bp8.2).0.5/Kme; 


end %end I loop 
end %end J loop 


% Table inputs for average actuator gains vs inclination follows: 
Inclination=[0.01 10 20 30 40 50 60 70 80 90 100 110 120]*pi/180; 
gx=[.967 .955 .922 .876 .826 .78 .739 .709 .691 .686 .691 .709 .739]; 
gy=[0.0001 .0678 .256 .46 .632 .762 .857 .923 .965 .981 .965 .923 .857]; 
gz=[.804 .781 .711 .614 522 .446 .39 .353 335 .333 .335 .353 .39]; 


toc 
NPS (NPSATACSINTRVW.MDL) 


See Figure A2. 
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Figure A2. NPSATACSIntRvw.mdl. 
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NPS (GRABDATASET.M) 


“%grabdataset.m written by Lt Todd Zirkle for manipulating data captured from 
% a SIMULINK run during analysis for NPSSAT1 NPS Monterey 
% 

% This program measures and plots the instantaneous and cumulative total of 
% commanded rod torques and the time required to achieve a 3-axis steady 

% state pointing accuracy of less than 10 degrees. 

% This program also solves and plots the area under the instantaneous rod 

% total curve to produce an energy index for comparison with other runs 

% using different initial conditions. 

% 

% Eleven variables are recorded in a matrix as a ‘.mat file’ from 

%  asimulation using one of two modified SIMULINK simulations originally 
% written by Prof Barry Leonard. 

%  ‘Lab17.mdl’ 05/18/01 and ‘NPSATACSIntRvw.mdl’ 07/16/01 


% This program picks off some of the desired variables. 

% The desired torque rod outputs in Am“2 are recorded as 

% ‘rodl’, ‘rod2’, and ‘rod3’ 

% The total commanded torque from the three rods combined is ‘rodtot’ 


% After loading the desired ‘.mat’ file containing the data matrix 
% into the MATLAB workspace, the matrix size is used to edit this 
%  program’s array sizes. 


clear ans xacksis DataArray rodone rod1 rod2 rod3 rodtot 
clear rodtwo rodthree rodtotal 

clear PhdT ThdT PsdT PhT ThT PsT 

clear Phd Thd Psd Ph Th Ps lasttime 

clear Timegap RulerTime miniarea areatot areaCumPlot 


DataArray = NPn7600'; % matrix name from workspace containing data 

xacksis = DataArray(1:11370,1); % array size must be updated to 
% match matrix size 

xaxiss=xacksis'; 

rod1= DataArray(1:11370,2); rodone = rod1'; 

rod2= DataArray(1:11370,3); rodtwo = rod2'; 

rod3= DataArray(1:11370,4); rodthree = rod3’; 

rodtot= DataArray(1:11370,5); rodtotal = rodtot’; 

PhdT= DataArray(1:11370,6); Phd = PhdT'; 

ThdT= DataArray(1:11370,7); Thd = ThdT’; 

PsdT= DataArray(1:11370,8); Psd = PsdT'; 

PhT= DataArray(1:11370,9); Ph = PhT'; 

ThT= DataArray(1:11370,10); Th = ThT’; 

PsT= DataArray(1:11370,11); Ps = PsT’; 
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areatot=0; % essentially starts the areatot at the origin which is zero 


for k = 1:11369 % will be one less than the size of the matrix dataset 
Timegap(k)= xacksis(k+1)-xacksis(k); % delta t for integrating area under curve 
holdtotal=areatot; 
miniarea(k)=rodtot(k)*Timegap(k); % represents area under curve 
% from (t+n) to (tt+n+timegap) 
areatot=holdtotal+miniarea(k); % cumulative area up til this point 
areaCumPlot(k)=areatot; % records cumulative area at each time interval 
end 


areaCumPlot(11370)=areaCumPlot(11369); % the last delta t is reused for 
% the last miniarea calculation 
lasttime=xacksis(11370); % last time of simulation recorded for plot label use 


figure % plots rod total commanded torques vs. time 
plot(xaxiss,rodtot,'r’),title(NPn (Lab17) Rod Total vs. Time’), 
ylabel('(Am“2)'), grid 
xlabel('Time (sec)') 


figure % plots cum area under curve vs. time [energy index (Am‘2 * sec)] 
plot(xaxiss,areaCumPlot,'r'), 

title({'NPn(Lab17) area Rod Total ',;num2str(areatot),' in ';num2str(lasttime),' sec']) 
ylabel(''), grid 

xlabel('Time (sec)') 


NPS (CREATE_VARIABLES.M) 


“%ocreate variables — written by Lt Todd Zirkle for capturing data from 

% NRL program written by Dr. Glenn Creamer, NRL Code 8231 
% (nps.m) 

% This program captures 11 variables for analysis and plotting purposes 

% Before running this program, run ‘nps.m’ after ensuring the initial conditions 
% have been set to the desired values. 

% This program will capture the variables by assigning them 

% to new variable names. 

% After renaming the 11 variables, the program clears all ‘nps.m’ variables. 

% The resulting data is saved by saving the MATLAB workspace 

% toa data folder. 

% It is recommended that the ‘.mat’ file’s name in some way reflects 

% the initial conditions used to produce the data, since the only other 

% way to ascertain the initial conditions is to pick out the first row of 

% data from the six variables assigned to rates and angles. 

% EXAMPLE FILE NAME: ‘bdot_pNp_testcase1.mat’ 

% In this example, pNp stands for 

% ‘phidot’ =.01 rad/s, ‘thdot’ =-.1 rad/s, ‘psidot’ =.01 rad/s 
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clear brod1 brod2 brod3 brodtot bww1 bww2 bww3 bpitt brole byeaw 


brod1=[m1]; % command to torque rod 1 (phi) ((from 0 to 30 Am‘2)) 
brod2=[m2]; % command to torque rod 2 (theta) ((from 0 to 30 Am“2)) 
brod3=[m3]; % command to torque rod 3 (psi) ((from 0 to 30 Am‘2)) 
brodtot=[rod3tot]; % rod3tot = the sum of [abs(rod1)+abs(rod2)+abs(rod3) | 

bwwl=[w1]; 

bww2=[w2]; 

bww3=[w3]; 


bpitt=[pitch_angle]; 
brole=[roll_angle]; 
byeaw=[yaw_ angle]; 
btime=[time]; 


clear yaw xyz_from_orb xyz_from_eciu_ext u_gg u_mag 

clear watt_hrs total t t_final theta true_anomaly 

clear store_count state sqrt2 satellite inertia inv_satellite_inertia 
clear ra radius rates roll r_ecir_pqwr_xyzr_xyz_unit 

clear pitch quats r2d r_ecef orb_from_eci orb_from_pqw 

clear orbit_rate period inv_satellite k k_ mag m_max 

clear eci_from_pqw greenwich_time i inc earth_gravity_constant 
clear earth radius earth_rate ecef_from eci d2r dipole 

clear dipole_off dt count bdot bfield_ecef bfield_xyz 

clear bfield_xyz_old ans ampsec_total altitude 

clear ecef_from_eci 

clear time wl w2 w3 w4 ml m2 m3 b1 b2 b3 

clear rate_goal rg time_counter pwrl pwr2 pwr3 pwr4 pwrSS 
clear pwr6 pwr7 pwr8 percent _pwr rod_current rod3tot 

clear pwr_u ugg! ugg2 ugg3 umagl umag2 umag3 

clear yaw_angle pitch_angle roll_angle libration 











NPS (NPSSATITESTCASE.M) 


%NPSsat1 Test case data manipulator 

% 

% Written by Lt Todd Zirkle for manipulating data captured from 

% MATLAB programs and SIMULINK runs during analysis 

% for NPSSAT1 NPS Monterey 

% 

% This program computes and plots the results of combining 

% two different simulations’ data. A Bdot control law simulation is 

% used to arrest rotation rates up to a predetermined rate goal, then the 
% final states of that stimulation are used for the initial conditions of the 
% second program, which further arrests the rotation rates and then 

% achieves a predetermined three-axis pointing accuracy. 
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% 
% 
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Uses an NRL program written by Dr. Glenn Creamer, NRL (nps.m) 
Uses two modified SIMULINK simulations originally written 
by Prof Barry Leonard. ‘Lab17.mdl’ 05/18/01 
‘NPSATACSIntRvw.mdl’ 07/16/01 


This program essentially combines two other programs written by Lt Zirkle: 
‘grabdataset.m’ and ‘create_variables.m’ 
As the original programs gathered rodtotal and energy index vs. time plots 
for the Dr. Creamer program and the Prof Leonard programs, this program 
does the same for the case of both programs being run in a handoff fashion 
as described above. 


This program captures torque rod total and time variables for analysis and 
plotting purposes 

Before running this program, run ‘nps.m’ after ensuring the initial conditions 
have been set to the desired values. 

This program will capture the variables by assigning them to new variable 
names. After renaming the variables, the program clears all non-essential 
variables. 


This program measures and plots the instantaneous and cumulative total of 
commanded rod torques and the time required to achieve a 3-axis steady 
state pointing accuracy of less than 10 degrees. 

This program also solves and plots the area under the instantaneous rod total 
curve to produce an energy index for comparison with other runs using 
different initial conditions. 


After loading the desired ‘.mat’ files containing data matrices into the 
MATLAB workspace, the matrix sizes are used to edit this program’s array 
sizes. 

The data are saved as a ‘.mat’ file. One of the two Prof Leonard programs 
is run after initially running the Bdot program. This data is stored as 
“handoff data” 


first import handoff data, then bdot data 
ensure handoff variable name appears on the 
line that creates the HDataArray matrix 
then ensure the matrix sizes are updated 


clear Hxacksis HDataArray Hrodone Hrod1 Hrod2 Hrod3 Hrodtot 
clear Hrodtwo Hrodthree Hrodtotal 

clear HPhdT HThdT HPsdT HPhT HThT HPsT 

clear HPhd HThd HPsd HPh HTh HPs 

clear HTimegap HRulerTime Hminiarea Hareatot HareaCumPlot 
clear newtime Composite_Rodtotal Composite_Time bdottime 
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HDataArray = NPnhandoff6600'; % input name of handoff matrix 
Hxaxis = HDataArray(1:8135,1); % input size of handoff matrix 

% Hxaxis will be the (time (x axis)) 
Hxacksis = Hxaxis'; 
bdottime = btime'’; % time (x-axis) before handoff 
Hrodtot= HDataArray(1:8135,5); Hrodtotal = Hrodtot'’; 
HTimegap(1)= Hxacksis(2)-Hxacksis(1); 
newtime(1)=6908; % final bdot time must be entered here 


for k=1:1151 % input matrix size of bdot matrix 
Composite Rodtotal(k)=brodtot(k); 
Composite_Time(k)=bdottime(k); 

end 


Composite Rodtotal(1152)=Hrodtot(1); 
Composite_Time(1152)=newtime(1); 


for k = 2:8134 % will be one less than handoff matrix size 
HTimegap(k)= Hxacksis(k+1)-Hxacksis(k); 
newtime(k)=newtime(k-1)+HTimegap(k); 
Composite _Rodtotal(1151+k)=Hrodtot(k); % input one less than bdot size 
Composite_Time(1151+k)=newtime(k); —_% input one less than bdot size 
end 





HTimegap(8135)= Hxacksis(8135)-Hxacksis(8134); % h-off x2, h-off less one 

newtime(8135)=newtime(8134)+HTimegap(8135); % h-off, h-off less one, h-off 
Composite _Rodtotal(1151+8135)=Hrodtot(8135); | % bdot plus h-off size, h-off 
Composite_Time(1151+8135)=newtime(8135); % bdot plus h-off size, h-off 


figure % plots composite rod total vs. time 
plot(Comp_Time,Comp_Rodtotal,'r’),title("Handoff Rod Total vs Time’), 
ylabel(' Am“2 '), grid 

xlabel('Time (sec)') 


NPS (RODTOTAL_VS_TIME.M) 


%Rodtotal vs Time written by Lt Todd Zirkle to plot simulated torque rod data 
% for analysis of simulations supporting NPSSAT1 NPS Monterey 


% This program captures a row of data from a saved 11 row matrix 

% The matrix was produced from a simulation using one of two modified 
% SIMULINK simulations originally written by Prof Barry Leonard. 

% ‘Lab17.mdl’ 05/18/01 and ‘NPSATACSIntRvw.mdl’ 07/16/01 
% This program could be used to capture any row for plotting purposes 
% This program grabs torque rod totals and plots them vs. time 
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clear ans HDataArray Hxaxis Hxacksis Hrodtot 


HDataArray = NPnCombo6908'; % input matrix name from ‘.mat’ file 

% loaded into workspace prior to run 
Hxaxis = HDataArray(1:1720,1); % input array size here from matrix size 
Hxacksis = Hxaxis'; % time axis 
Hrodtot= HDataArray(1:1720,5); Hrodtotal = Hrodtot'’; 


figure % plots torque rod total (sum of absolute values of 3 rods) vs time 
plot(Hxaxis,Hrodtot,'r'’),title((B.Leonard Combo with Bdot Rod Total vs Time’), 
ylabel(' Am*2 '), grid 


xlabel('Time (sec)') 
NPS (PLOTRESULT2.M) 
“%plotresult2.m written by Todd A. Zirkle for graphical analysis 
% of NPSSATI1 data) NPS Monterey 
figure 


subplot(411),plot(time,w1), 

title(['w_t_o t_a 1<'num2str(rate_goal),' degrees ('snum2str(rg),' rad) in 
"num2str(time_counter),' sec']), 

ylabel('‘w_x(rad/sec)'),grid 
subplot(412),plot(time,w2),ylabel('w_y(rad/sec)'),ylabel('w_ y(rad/sec)'), grid 
subplot(413),plot(time,w3),ylabel('w_z(rad/sec)'),ylabel('w_z(rad/sec)'), grid 
subplot(414),plot(time,w4),ylabel(‘w_t_o t_a_1'),ylabel(‘w_t_o t_a_1'),grid 
xlabel(‘time(sec)') 


figure 

subplot(511),plot(time,pwr_u),title(['Total Power =',num2str(pwr3),' 
W-hrs']),ylabel(' W-hrs'), grid 
subplot(512),plot(time,watt_hrs_total),ylabel(' W-hrs'),ylabel(' W-hrs'), grid 
subplot(513),plot(time,m1,'r'),ylabel(‘Am‘’2'), ylabel('‘Am“2'), grid 
subplot(514),plot(time,m2,'g'),ylabel('Am%2'), ylabel('‘Am“2'), grid 
subplot(515),plot(time,m3,'b'),ylabel('Am%2'), ylabel('‘Am”2'), grid 
xlabel(‘time(sec)') 


figure 

subplot(511),plot(time,pwrSS),title(['SS Power = 'snum2str(pwr6),' 
W-hrs']),ylabel(' W-hrs'), grid 
subplot(512),plot(time,pwr6),ylabel(' W-hrs'),ylabel(’ W-hrs'), grid 
subplot(513),plot(time,m1,'r'),ylabel('Am’2'), ylabel('‘Am“2'), grid 
subplot(514),plot(time,m2,'g'),ylabel('Am%2'), ylabel('‘Am”2'), grid 
subplot(515),plot(time,m3,'b'),ylabel('Am%2'), ylabel('‘Am“2'), grid 
xlabel('time(sec)') 
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figure 

subplot(411),plot(time,rod3tot),title({'3 Rod torque = 'snum2str(rod3tot),' Am*2']), 
ylabel(' Am‘’2'), grid 

subplot(412),plot(time,m1,'r'),ylabel('Am’2'), ylabel('‘Am“2'), grid 
subplot(413),plot(time,m2,'g'),ylabel('Am%2'), ylabel('‘Am”2'), grid 
subplot(414),plot(time,m3,'b'),ylabel('Am%2'), ylabel('‘Am“2'), grid 
xlabel('time(sec)') 


figure 


plot(time,rod3tot),title({'3 Rod torque = 'snum2str(rod3tot),' Am“%2']), 
ylabel(' Am“2'), grid 


90 


BIBLIOGRAPHY 


Franklin, G. F., Powell, J. D., Emami-Naeini, A., “Feedback Control Of Dynamic 
Systems.” Addison Wesley, 1988. 


Martel, F., Pal, P. K., and Psiaki, M., “Active Magnetic Control System for 
Gravity Gradient Stabilized Spacecraft.” Second Annual AIAA/USU Conference on 
Small Satellites, (26-28 September 1988, Logan, UT). Washington, DC: AIAA. pg 7, 13 

“MATLAB Control System User’s Guide.” The Math Works, Natick, MA. 1996. 

Psiaki, M. L., Martel, F., and Pal, K. P., ““Three- Axis Attitude Determination via 
Kalman Filtering of Magnetometer Data”, Journal of Guidance, Control and Dynamics, 


Vol. 23, No. 3, 1990.pp 506-514. 


Riley, F. E., Sailor, J. D., “Space Systems Engineering”, McGraw-Hill Book 
Company, Inc., 1962.pp 133-151. 


Sidi, M. J., “Spacecraft Dynamics and Control.” Cambridge University Press. 
1997. 


Wertz, J. R. “Spacecraft Attitude Determination and Control.” Kluwer Academic 
Publishers. Boston, MA. 1994. 


Wertz, J. R., Larson, W. J., “Space Mission Analysis and Design.” 3rd Edition, 
Microcosm Press. 1999. 


91 


THIS PAGE INTENTIONALLY LEFT BLANK 


92 


INITIAL DISTRIBUTION LIST 


Defense Technical Information Center ........... cece ccc eccccccccccccccceccccceceeceseecececeeeeses 
8725 John J. Kingman Road, Suite 0944 
Ft. Belvoir, VA 22060-6218 


Dudley hex WiDrary wc iceys oasodeedian neste idee MGs aula ak aha 
Naval Postgraduate School 

411 Dyer Road 

Monterey, CA 93943-5101 


93 


